normalSample
Statistics Mar10 Timm Revised

How to Sample from a Normal Distribution]

Problem
=======

Given the mean "_m_" and standard deviation "_s_" of a normal distribution,
you want to generate "_n"_ numbers from that distribution.

Solution
========

Using the _box\_muller_ transform. That is,
take two samples from the uniform distribution on the interval
(0, 1\] and map them to two normally distributed samples. The polar
form takes two samples from a different interval, \[-1, +1\], and maps
them to two normally distributed samples without the use of sine or
cosine functions.

No, I don't really understand it either. 

 function nsample(m,s,n,a,     i) {
     for(i=1;i<=n;i++)
         a[i]=normal(m,s)
 }
 function normal(m,s) {
     return m+box_muller()*s;
 }
 function box_muller(m,s,    n,x1,x2,w) {
     w=1;
     while (w >= 1) {
         x1= 2.0 * rand() - 1;
         x2= 2.0 * rand() - 1;
         w = x1*x1 + x2*x2};
     w = sqrt((-2.0 * log(w))/w);
     return x1 * w;
 }

Example
=======

The following code should report means and standard
deviations near 100 and 10 (respectively).

 function Nsample(   a,i,m,samples)  {
     seed(1)
     samples=1000
     nsample(100,10,samples,a)
     for(i in a)
         nkeep(a[i],m)
     print nmean(m), nsd(m);
 }

Author
======

Tim Menzies
